/* This program calculates the sup-t critical value for a set of coefficients
in a regression, then posts the critical value to a set of stored results.
This is designed to be called by estadd from the estout package*/

program define st_crit, rclass

syntax varlist(fv ts) [,  n(integer 10000) alpha(real 95) ]

tempname oldresults b V b0 crit




/* Pull coefficient vector and VCE matrix*/
mat `b' = e(b)
mat `V' = e(V)
local newnames `:colnames `b''

/* Check that the given varlist are all in the results*/
foreach var of varlist `varlist'{
	if regexm("`: colnames `b''","`var'") ==0 {
		display as error "`var' not found in estimation results*"
		error 198
	}
	/* Assuming they are, append an equation name to make extraction easy*/
	else {
		local newnames = regexr("`newnames'","((^`var')|( `var'))"," extract:`var'") //To avoid duplicates, either require match at the beginning, or a space before
	}
}
matrix colnames `V' = `newnames'
matrix rownames `V' = `newnames'

/* Extract just rows and columns from V for the variables of interest,
and convert to a correlation matrix (so it's unit-independent)*/


mat `V' = corr(`V'["extract:","extract:"])

/* Max a list of tempvar names for the random draws*/
local newlist
forvalues i = 1/`: word count `varlist''{
	tempvar draw`i'
	local newlist `newlist' `draw`i''
}

/* Make N draws from a multivariate normal for each variable, using the betas and VCE we extracted*/
preserve

mat `b0' = J(1,`: word count `varlist'',0)
qui drawnorm `newlist', means(`b0') corr(`V') n(`n') double clear

/* Now extract the max of the absolute values of the betas from each draw*/
foreach var of varlist `newlist'{
	qui replace `var' = abs(`var') 
}

tempvar maxdraw
qui egen `maxdraw' = rowmax(`newlist')

/*Critical value is then the relevant percentile of the distribution of max's*/

_pctile `maxdraw' , percentiles(`alpha')


return scalar crit = [100 * (1-2*(1-normal(r(r1))))]

end
